   
# ---------------------------------------------------------------
## Scripts to load CITUP and PyClone to plot trees (Figure 2)
## Author: Veronique LeBlanc
## https://github.com/vleblanc/GBM-PDO-paper/blob/main/WES/fig2_trees.R
## https://github.com/hanhanial/PyClone_VI/blob/c2b7af0eb01d2f1e5fcadf907948bec3b5cb98d1/scripts/S05_plotting_citup-clone-trees.r
# ---------------------------------------------------------------

#library(ggplot2)
library(argparser)
library(reshape2)
library(data.table)
library(dplyr)
library(parallel)

##############################################################################

argp <- arg_parser("Plot the CitupPlot")
argp <- add_argument(argp, "--smg_gene_file" , help="")
argp <- add_argument(argp, "--output_path" , help="")
argp <- add_argument(argp, "--sample_file" , help="")
argp <- add_argument(argp, "--class_order" , help="")
argp <- add_argument(argp, "--class_a_order" , help="")
argp <- add_argument(argp, "--clone_t" , help="")
argp <- add_argument(argp, "--ccf_file" , help="")

argv <- parse_args(argp)

smg_gene_file <- argv$smg_gene_file
output_path <- argv$output_path
class_order <- argv$class_order
class_a_order <- argv$class_a_order
sample_file <- argv$sample_file
clone_t <- as.numeric(argv$clone_t)
ccf_file <- argv$ccf_file

if(1!=1){

  work_dir <- "~/20220915_gastric_multiple/dna_combinePublic"
  smg_gene_file <- paste0( work_dir , "/mutsig_check/report_chooseclone.list")
	output_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/selectGCClone"
  ccf_file <- paste(work_dir,"/mutationTime/result/All_CCF_mutTime.addShare.tsv",sep="")
  clone_t <- 0.6
  class_order <- paste(work_dir,"/config/Class_order.list",sep="")
  class_a_order <- paste(work_dir,"/config/Class_order_sub.list",sep="")
  sample_file <- paste(work_dir,"/config/tumor_normal.class.list",sep="")

}

dir.create(output_path , recursive = T)

##############################################################################
# 已报道驱动基因
dat_smg <- data.frame(fread(smg_gene_file , header = F))
colnames(dat_smg) <- "Gene_Symbol"

##############################################################################
class_order <- data.frame(fread(class_order))
class_a_order <- data.frame(fread(class_a_order))
info <- data.frame(fread(sample_file))
dat_ccf <- data.frame(fread(ccf_file))

## 5个样本只用于最后观察
#info <- subset( info , Type != "IM + IGC + DGC" )

##############################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
dat <- subset( dat_ccf , Variant_Classification %in% Variant_Type )

##############################################################################
## 判断每个基因多少突变为IM和IGC或DGC共享、IGC和DGC中主克隆突变的占比
## IGC和DGC的IM各自独有的突变数目

igc_sample <- unique(subset( info , Type=="IM + IGC"  )$ID)
dgc_sample <- unique(subset( info , Type=="IM + DGC" )$ID)
igc_dgc_sample <- unique(subset( info , Type=="IM + IGC + DGC" )$ID)

result <- Reduce(function(x,y)bind_rows(x,y),mclapply( unique(dat$Hugo_Symbol) , function(geneN){

  tmp_ccf <- subset( dat , Hugo_Symbol == geneN )
  tmp_ccf_igc <- subset( tmp_ccf , ID %in% igc_sample )
  tmp_ccf_dgc <- subset( tmp_ccf , ID %in% dgc_sample )
  tmp_ccf_igc_dgc <- subset( tmp_ccf , ID %in% igc_dgc_sample )

  ## 判断每个基因多少突变为IM和IGC或DGC共享
  shareNum_igc <- length(unique(subset( tmp_ccf_igc , Share=="TRUE" )$ID))
  shareNum_dgc <- length(unique(subset( tmp_ccf_dgc , Share=="TRUE" )$ID))

  shareNum_igc_igcdgc <- length(unique(subset( tmp_ccf_igc_dgc , Share=="TRUE" & Class %in% c("IGC") )$ID))
  shareNum_dgc_igcdgc <- length(unique(subset( tmp_ccf_igc_dgc , Share=="TRUE" & Class %in% c("DGC") )$ID))

  ## IGC和DGC的IM多少突变为独有
  uniqueIMNum_igc <- length(unique(subset( rbind( tmp_ccf_igc , tmp_ccf_igc_dgc) , Share=="FALSE" & Class == "IM" )$ID))
  uniqueIMNum_dgc <- length(unique(subset( rbind( tmp_ccf_dgc , tmp_ccf_igc_dgc) , Share=="FALSE" & Class == "IM" )$ID))

  ## IGC和DGC中主克隆突变的占比
  tmp_ccf_igc_unique <- subset( rbind( tmp_ccf_igc , tmp_ccf_igc_dgc) , Class == "IGC" )
  tmp_ccf_igc_unique <- tmp_ccf_igc_unique %>%
  group_by( ID ) %>%
  summarize( CCF_adj = max(CCF_adj) )
  igc_im_cloneRatio <- length(unique(subset( tmp_ccf_igc_unique , CCF_adj >= clone_t )$ID))/length(unique(tmp_ccf_igc_unique$ID))

  tmp_ccf_dgc_unique <- subset( rbind( tmp_ccf_dgc , tmp_ccf_igc_dgc) , Class == "DGC" )
  tmp_ccf_dgc_unique <- tmp_ccf_dgc_unique %>%
  group_by( ID ) %>%
  summarize( CCF_adj = max(CCF_adj) )
  dgc_im_cloneRatio <- length(unique(subset( tmp_ccf_dgc_unique , CCF_adj >= clone_t )$ID))/length(unique(tmp_ccf_dgc_unique$ID))

  ## GC中的克隆占比
  tmp_ccf_unique <- subset( tmp_ccf , Class %in% c("IGC" , "DGC") )
  tmp_ccf_unique <- tmp_ccf_unique %>%
  group_by( ID ) %>%
  summarize( CCF_adj = max(CCF_adj) )
  im_cloneRatio <- length(unique(subset( tmp_ccf_unique , CCF_adj >= clone_t )$ID))/length(unique(tmp_ccf_unique$ID))

  ## 计算该基因在多少IGC/DGC中不与已知驱动基因共同出现
  mut_igc_sample <- c(unique(tmp_ccf_igc$ID) , unique(tmp_ccf_igc_dgc$ID))
  mut_dgc_sample <- c(unique(tmp_ccf_dgc$ID) , unique(tmp_ccf_igc_dgc$ID))

  # 计算在多少样本该基因与共享驱动基因一起出现
  tmp <- subset( dat , ID %in% mut_igc_sample & Hugo_Symbol %in% dat_smg$Gene_Symbol & Class == "IGC")
  tmp <- tmp %>%
  group_by( ID ) %>%
  summarize( Hugo_Symbol = paste0(Hugo_Symbol , collapse = ",") )
  unique_num <- length(mut_igc_sample) - length(unique(tmp$ID)) 
  unique_igc_num <- unique_num

  tmp <- subset( dat , ID %in% mut_dgc_sample & Hugo_Symbol %in% dat_smg$Gene_Symbol & Class == "DGC")
  tmp <- tmp %>%
  group_by( ID ) %>%
  summarize( Hugo_Symbol = paste0(Hugo_Symbol , collapse = ",") )
  unique_num <- length(mut_dgc_sample) - length(unique(tmp$ID)) 
  unique_dgc_num <- unique_num

  tmp <- subset( dat , ID %in% c(mut_dgc_sample , mut_igc_sample) & Hugo_Symbol %in% dat_smg$Gene_Symbol & Class %in% c("IGC" , "DGC") )
  tmp <- tmp %>%
  group_by( ID ) %>%
  summarize( Hugo_Symbol = paste0(Hugo_Symbol , collapse = ",") )
  unique_num <- length(unique(c(mut_dgc_sample , mut_igc_sample))) - length(unique(tmp$ID)) 
  unique_gc_num <- unique_num

  tmp_res <- data.frame( Hugo_Symbol = geneN , 
    shareNum_igc = shareNum_igc + shareNum_igc_igcdgc , shareNum_dgc = shareNum_dgc + shareNum_dgc_igcdgc , shareNum_gc = shareNum_dgc + shareNum_dgc_igcdgc  + shareNum_igc , # 多少样本属于IM和GC共享
    uniqueIMNum_igc = uniqueIMNum_igc , uniqueIMNum_dgc = uniqueIMNum_dgc , uniqueIMNum_gc = uniqueIMNum_igc + uniqueIMNum_dgc , # 多少样本属于IM独有
    igc_cloneRatio = igc_im_cloneRatio , dgc_cloneRatio = dgc_im_cloneRatio , gc_cloneRatio = im_cloneRatio ,# 胃癌中多少比例为clone
    NoDriverConoccur_igc = unique_igc_num , NoDriverConoccur_dgc = unique_dgc_num , NoDriverConoccur_gc = unique_gc_num
    )

  tmp_res

},mc.cores=20))

##############################################################################
## 筛选标准
## 2、IM中不存在，GC中该基因超过80%均为主克隆：

## 已报道的驱动
gene_list_report <- unique(subset(result , Hugo_Symbol %in% dat_smg$Gene_Symbol )$Hugo_Symbol)

## 1、IM中存在，GC中必定保留，且GC中该基因超过80%均为主克隆
gene_list1 <- unique(subset(result , shareNum_gc >= 2 & uniqueIMNum_gc == 0 & gc_cloneRatio >= 0.8 & NoDriverConoccur_gc >= 2)$Hugo_Symbol)

## 2、IM中不存在，GC中必定保留，且GC中该基因均为主克隆
gene_list2 <- unique(subset(result , shareNum_gc == 0 & uniqueIMNum_gc == 0 & gc_cloneRatio == 1 & NoDriverConoccur_gc >= 2)$Hugo_Symbol)

## 3、IM存在GC未保留的突变，IGC/DGC特异型驱动突变
# 仅在IGC或DGC的IM中独立存在，在另一种胃癌亚型中均为主克隆突变
gene_list3_igc <- unique(subset(result , shareNum_igc > 0 & uniqueIMNum_igc == 0 & shareNum_dgc == 0 & uniqueIMNum_dgc > 0 & igc_cloneRatio == 1 & NoDriverConoccur_igc >= 2)$Hugo_Symbol)
gene_list3_dgc <- unique(subset(result , shareNum_dgc > 0 & uniqueIMNum_dgc == 0 & shareNum_igc == 0 & uniqueIMNum_igc > 0 & dgc_cloneRatio == 1 & NoDriverConoccur_dgc >= 2)$Hugo_Symbol)

result2 <- result
result2$CloneChoose <- "NoChoose"
result2$CloneChoose <- ifelse( result2$Hugo_Symbol %in% gene_list_report , "ReportDriver" , result2$CloneChoose )
result2$CloneChoose <- ifelse( result2$Hugo_Symbol %in% gene_list1 , "IM_Share" , result2$CloneChoose )
result2$CloneChoose <- ifelse( result2$Hugo_Symbol %in% gene_list2 , "GC_Private" , result2$CloneChoose )
result2$CloneChoose <- ifelse( result2$Hugo_Symbol %in% gene_list3_dgc , "DGC_Choose" , result2$CloneChoose )
result2$CloneChoose <- ifelse( result2$Hugo_Symbol %in% gene_list3_igc , "IGC_Choose" , result2$CloneChoose )

images_name <- paste0(output_path , "/GCClone_gene.reord.tsv" )
write.table(result2 , images_name , sep = '\t' , quote=F , row.names=F)
